
#include "slepdef.h"
#include "slepblas.h"
#include <string.h>
#include <math.h>
#include <assert.h>
#include <malloc.h>


void slep_daAxpby(slep_transposity trans, double alpha,
				  slep_matrix* A, slep_vector* x,
				  double beta, slep_vector* y){
	static int i,j;
	static double tempd=0;

	if(alpha == 0){
		slep_dscal(beta,y);
		return;
	}

	if(A->sparsity == SlepDense && x->sparsity == SlepDense && y->sparsity == SlepDense){
		if( trans == SlepNoTranspose){
			assert( A->ld >= A->m );
			assert( A->n == x->dim );
			assert( A->m == y->dim );
			for(i=0; i < A->m ; i++){
				y->data.vctr[i] *= beta;
				for(j=0; j < A->n; j++)
					y->data.vctr[i] += alpha * A->data.mtrx[i + j*(A->ld)] * x->data.vctr[j] ;
			}
		}
		else{
			assert( A->ld >= A->m );
			assert( A->n == y->dim );
			assert( A->m == x->dim );
			for(i=0; i < A->n ; i++){
				y->data.vctr[i] *= beta;
				for(j=0; j < A->m; j++)
					y->data.vctr[i] += alpha * A->data.mtrx[i*(A->ld) + j] * x->data.vctr[j] ;
			}
		}
		return;
	}

	if( A->sparsity==SlepSparse && x->sparsity==SlepDense && y->sparsity==SlepDense){

		/* This wont be efficient TODO conjoin with the code below */
		slep_dscal(beta,y);

		if ((A->data.spmtrx.format == SlepCompressedSparseColumn) 
			&& (trans == SlepNoTranspose) )
			for( i=0 ; i < A->n ; ++i )
				for( j= A->data.spmtrx.Index[i]; j < A->data.spmtrx.Index[i+1];j++)
					y->data.vctr[(A->data.spmtrx.ColumnsORRows[j])] += (alpha * (A->data.spmtrx.values[j])*(x->data.vctr[i]) ); 
		else if ((A->data.spmtrx.format == SlepCompressedSparseRow) 
			&& (trans == SlepTranspose) )
			for( i=0 ; i < A->m ; ++i )
				for( j= A->data.spmtrx.Index[i]; j < A->data.spmtrx.Index[i+1];j++)
					y->data.vctr[(A->data.spmtrx.ColumnsORRows[j])] += (	alpha * (A->data.spmtrx.values[j])*(x->data.vctr[i]) );
		else if ((A->data.spmtrx.format == SlepCompressedSparseRow) 
			&& (trans == SlepNoTranspose) )
			for(i=0; i < A->m; ++i)
				for(j= A->data.spmtrx.Index[i]; j < A->data.spmtrx.Index[i+1];j++)
					y->data.vctr[i] += alpha * (A->data.spmtrx.values[j]) * (x->data.vctr[A->data.spmtrx.ColumnsORRows[j]]) ;
		else if ((A->data.spmtrx.format == SlepCompressedSparseColumn) 
			&& (trans == SlepTranspose) )
			for(i=0; i < A->n; ++i)
				for(j= A->data.spmtrx.Index[i]; j < A->data.spmtrx.Index[i+1];j++)
					y->data.vctr[i] += alpha * (A->data.spmtrx.values[j]) * (x->data.vctr[A->data.spmtrx.ColumnsORRows[j]]) ;
		return;
	}

	if(A->sparsity == SlepDense && x->sparsity == SlepSparse && y->sparsity == SlepDense){
		if( trans == SlepNoTranspose){
			assert( A->m  == y->dim );
			assert( x->data.spvctr.index[x->data.spvctr.nnz - 1] <= A->n);
			for(i=0;i< A->m ; i++){
				y->data.vctr[i* y->inc] *= beta;
				for(j=0, tempd=0; j< x->data.spvctr.nnz; j++)
					tempd += (A->data.mtrx[(x->data.spvctr.index[j]) * (A->ld) + i]) * x->data.spvctr.values[j] ;
				y->data.vctr[i* (y->inc)] += alpha *tempd;
			}
		}
		else if( trans == SlepTranspose){
			assert( A->n  == y->dim );
			assert( x->data.spvctr.index[x->data.spvctr.nnz - 1] <= A->m);
			for(i=0;i< A->n ; i++){
				y->data.vctr[i* y->inc] *= beta;
				for(j=0, tempd=0; j< x->data.spvctr.nnz; j++)
					tempd += (A->data.mtrx[i * (A->ld) + (x->data.spvctr.index[j])]) * x->data.spvctr.values[j] ;
				y->data.vctr[i* (y->inc)] += alpha * tempd;
			}
		}
		return ;
	}

}

void slep_dtrmv(slep_transposity trans, slep_matrix* A,
				slep_vector* x){
	static int i,j;
	if(trans == SlepNoTranspose){
		if(A->uplobanded == SlepUpper)
			for(i=0; i < A->m; i++){
				x->data.vctr[i] *= A->data.mtrx[ (i*(i+3))/2 ];
				for(j=(i+1) ; j < A->m ; j++)
					x->data.vctr[i] += A->data.mtrx[((j*(j+1))/2) + i] * x->data.vctr[j];
			}
		else if(A->uplobanded == SlepLower)
			for(i= (A->m) - 1; i >= 0 ; i-- ){
				x->data.vctr[i] *= A->data.mtrx[ i + i*(A->m -((i+1)/2))] ;
				for(j=0; j < i ; j++)
					x->data.vctr[i] += x->data.vctr[j] * A->data.mtrx[i + j*(A->m -((j+1)/2))] ;
			}
/*		else if(A->uplobanded == SlepBandedUpper)
			;
		else if(A->uplobanded == SlepBandedLower)
			;*/
	}
	else if(trans == SlepTranspose){
		if(A->uplobanded == SlepUpper)
			for(j = (-1+ A->m); j >= 0 ; j--){
				x->data.vctr[j] *= A->data.mtrx[ (j*(j+3))/2 ];
				for(i= 0 ; i < j ; i++)
					x->data.vctr[j] += A->data.mtrx[((j*(j+1))/2) + i] * x->data.vctr[i];
			}
		else if(A->uplobanded == SlepLower)
			for(j= 0; j < A->m ; j++ ){
				x->data.vctr[j] *= A->data.mtrx[ j + j*(A->m -((j+1)/2))] ;
				for(i=j+1; i < A->m ; i++)
					x->data.vctr[j] += x->data.vctr[i] * A->data.mtrx[i + j*(A->m -((j+1)/2))] ;
			}
	}
	/* TODO Complete all the Banded part and also update for a full triangular matrix(current is for packed)*/
}

void slep_dtrsv(slep_transposity trans, slep_matrix* A,
				slep_vector* x){
}





